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The primary usefulness of the presented method is in the ability to 
represent the effects of high frequency linear response with accuracy, 
without requiring very small time steps in the analysis of dynamic 
response. The matrix exponential contains a series approximation to 
the dynamic model. However, unlike the usual modal analysis proce- 
dure which truncates the high frequency response, the approximation 
in the exponential matrix solution is in the time domain. By truncat- 
ing the series solution to the matrix exponential short, the solution 
is made inaccurate after a certain time. Yet, up to that time the 
solution is extremely accurate, including all high frequency effects. 
By taking finite time increments, the exponential matrix solution can 
compute the response very accurately. .Use of the exponential matrix 
in structural dynamics is demonstrated by simulating the free vibra- 
tion response of multi degree of freedom models of cantilever beams. 


INTRODUCTION 

The matrix exponential has been known in matrix theory as a method of solution for systems 
of differential equations. However, it has not been applied to the solution of structural dynamics 
problems. This method may be useful in some types of structural problems where modal decom- 
position is not practical or the number of modal vectors that can be accurately determined do not 
represent the true structural response. The matrix exponential contains a series approximation to 
the dynamic model. However, unlike the usual modal analysis procedure which truncates the high 
frequency response, the approximation in the exponential matrix solution is in the time domain. 
The exponential matrix method simulates the complete structural response, including the high 
frequency effects, but only for not too large values of the time parameter t . These properties make 
the exponential matrix method ideally suitable to complement the direct time-history integration 
of the equations of motion; to improve accuracy and to increase the integration stepsize. An 
updated Lagrangian formulation may be used at each integration step to recompute the matrix 
exponential with reference to the current state variables. 
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MATHEMATICAL BACKGROUND 

Before examining the implementation of the matrix exponential in the solution of structural 
dynamics equations, it is helpful to consider the homogeneous systems of differential equations. 

[X] = [A][X] (1) 

with the initial conditions 

[*(0)] = [I] ( 2 ) 

where [A] is an n x n constant matrix, [X] is an n x n matrix, the columns of which are individual 

unknown vectors and [I] is the n x n identity matrix. It can be shown that the fundamental 
solution [X] to Eq. (1) can be written as [1] 

oo 1 

[X } = e w t = j2h [A]ktk (3) 


where 


[A] fc = [A] ■ ■ • [A] (k terms) (4) 

It can also be shown that the infinite series given by Eq. (3) converges uniformly and absolutely 
for t in any bounded interval [1], Eq. (3) is referred to as the exponential matrix function. 

To s um marize the use of the matrix exponential in structural dynamics consider the equations 
of motion of a structure, written in the physical coordinates: 


[ M){x} + [C]{x} + [K}{x} = {F(t )} 


where the symbols have their usual meanings. By defining 
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the equations of motion, Eq. (5), can be writtten as a system of first order differential equations 
as 
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or in short notation redefining new symbols for the overall matrices and vectors in Eq. 
represent them as 


( 7 ) 

(7), we 


{*/} = [AJM + {/(<)} 

Introducing a new unknown vector {z(t)}, defined by the equation 

{y(t)} = [X(t)]{*(t)} 


( 8 ) 

( 9 ) 
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where [AT(f)] is the fundamental solution defined by Eqns. (1) and (3), we can write the time 
derivative of {y} as 

w = i*iw + m{i} (io) 

Substituting Eq. (10) into Eq. (8): 

[X}{z} + [X){z} = [A][X]{z} + {f} ( 11 ) 

Combining Eq. (1) with Eq. (11) we obtain 

[*]{*} = {/} ( 12 ) 

The solution of which can be written as [1] 

t 

{*(<)} = {*(0)} + J Wr)]- 1 {f(r)}dr (13) 

0 

and consequently 

t 

{»(*)} = WOIWO)} + [A-(*» /[jrwrq/o-)}*- (14) 

0 

with {^(0)} = {y(0)} since {t /} = [X]{z} and [X(0)] = [I], It can be shown that 

WOIPW = [X(t - t )] (15) 

if the matrix [.A] is independent of t. Thus 

t 

{»(<)} = W<)){y(0)} + J [X(t - T)]{f(r))dT (16) 

0 

where (y(^)} is the list of structural coordinate displacements and velocities as defined by Eq. (6) 

EXAMPLES 

Eq. (16) would give the correct solution for all t if the structure properties were independent 
of t and if [X(f)] could be computed with sufficient accuracy. In general, it is not practical to 
compute pf(<)] to a sufficient level of accuracy for Eq. (16) to be valid for all t. However, if t is 
not large, then Eq. (16) is expected to yield good results with relatively crude approximations 
[-^”(0] - •^' or example, Fig. 1 shows a comparison of the true free vibration response of a single 
degree of freedom system with the response computed by Eq. (16) when only 4 terms are included 
in Eq. (3). Eq. (16) matches the true response exactly, but only for approximately 1/2 period of 
vibration. Fig. 2 shows the same comparison when 16 terms are included in Eq. (3). In this case 
the exact simulation lasts for two periods of vibration. 
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Fig. 1.- Free Vibration Response of a Simple Oscillator 
(exponential series truncated after 4 terms) 



Fig. 2.- Free Vibration Response of a Simple Oscillator 
(exponential series truncated after 16 terms) 


Example with Two Degrees of Freedom 

As a two-degree-of-freedom physical example, consider a steel cantilever beam, two meters 
long, with a 0.1 meter square cross section, with discretized degrees of freedom as numbered in 
Fig. 3. 
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Fig. 3.- Example Beam with Coordinates 


The stiffness matrix for the numbered coordinates is assembled from beam element stiffnesses 
as: 


[K] = El 
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Or, we can represent the stiffness matrix [if] in terms of the submatrices as outlined by the dashed 
lines in Eq. (17) as 


[K] = 


K u 

K n 


I I<12 ' 

I K22 


(18) 


Using matrix condensation [2], the condensed stiffness matrix [7^*] corresponding to the first two 
coordinates is written as: 


UH = [Kn] - \K 17 }[K„)-'[K 21 ) 

or 


(19) 


[K 1 


El T 96 -30 

7 -30 12 


( 20 ) 


Substituting E = 200 GPa for steel and dropping the star from our notation for the condensed 
stiffness matrix, we obtain the stiffness matrix for the structure degree-of-freedom coordinates 
shown in Fig. 4 as 


[I<] = 


20 x 10 6 [ 96 
(7)(12) [-30 


-30 

12 


N/m 


( 21 ) 
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Fig. 4.- Coordinate numbering with condensed degrees-of-freedom 


Substituting p = 8, 000 kg/m 3 for the mass density of steel and considering only the transverse 
inertia of the beam, the lumped mass matrix for the two coordinates shown in Fig. 4 is written 
as: 


[M] 


80 0 
0 40 


kg 


( 22 ) 


Substituting these values for [M] and [K] in Eq. (5) and assuming that [C] and {/} are null, 
we define the undamped free vibration problem for this cantilever beam as: 


'80 

0 ' 

IM+io 6 ’ 
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40 

1 * 2 / L 


22.86 

-7.143 


-7.143] f*i\ f O'! 
2.857 ] \ x 2 J \ 0 J 


(23) 


The time history response due to a set of initial conditions can be obtained by a direct numer- 
ical integration procedure, such as the central difference method. In the following examination 
the central difference method is used to compare the performance of various exponential matrix 
solutions. 


Fig. 5 shows a comparison of the exact solution and a single step exponential matrix solution 
for the free vibration response of this example. In Fig. 5, only the first fifteen terms are included 
for the series solution of the matrix exponential as defined in Eq. (3). Fig. 6 shows the same 
comparison between the one step exponential matrix solution and the exact solution, but in this 
case thirty terms are used in the series definition of the matrix exponential. The solution depicted 
in Fig. 6 follows the exact solution for twice the time length as compared to the solution shown 
in Fig. 5, demonstrating the linear convergence of the exponential series. It is significant that 
the exponential matrix solution is identical to the true solution until very close to the divergence 
time. A second interesting point is that the simulation for both physical coordinates diverges 
simultaneously but in opposite directions. Fig. 7 shows the same comparison but considering 
twenty-nine terms in the series solution. When Fig. 7 is compared to Fig. 6, it is noted that 
divergence of the exponential simulation in these two figures go in opposite directions. It can be 
verified that divergence is always in a predictable direction, depending on having an odd or even 
number of terms in the series. The true benefit of the exponential matrix solution can be utilized 
when, for a given structural problem, the relationship between the accuracy of the exponential 
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Fig. 5 - Coordinate Displacements during Free Vibration. 
(Exponential Matrix Solution is Based on 15 Terms.) 




Time (Sec.) 

Fig. 6.- Coordinate Displacements during Free Vibration. 
(Exponential Matrix Solution is Based on 30 Terms.) 
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Fig. 7.- Coordinate Displacements during Free Vibration. 
(Exponential Matrix Solution is Based on 29 Terms.) 


matrix and the time duration of the accuracy of the solution can be predicted. Once the numerical 
technique for making this prediction is established, the exponential matrix solution can be used 
in discrete steps for extremely accurate time history analysis of dynamic systems. 

Next, it is informative to study the same problem, during a longer time interval. This time, 
exponential matrix solutions are combined, step by step, to render an extraordinarily accurate 
solution of the system. The step by step exponential matrix solution can be used to obtain an 
exact solution of a dynamic system, even if an extremely large stepsize is used. 

Figs. 8 and 9 show the exponential matrix solution and the central difference solution, 
respectively. In both figures the solid lines represent the “true” solution based on a central 
difference solution using a time increment of 0.0001 sec. The total simulation time is 0.1 sec. The 
exponential matrix solution is based on thirty terms in the series approximation. This exponential 
matrix solution with thirty terms diverges if the time interval is taken to be more than 0.021 sec., 
as depicted in Fig. 6. Accordingly, 0.020 sec is selected as the stepsize for the exponential 
matrix solution. The exponential matrix solution at each 0.02 sec time increment is plotted in 
Fig. 8 using square markers. The central difference solution diverges if the time interval is taken 
more than 0.0033 sec. In Fig. 9 the stepsize in the central difference solution is taken as 0.003 
sec. In comparing Figs. 8 and 9, the dramatic difference between the exponential matrix and 
the central difference step-by-step solutions illustrate the relative effectiveness of the exponential 
matrix method. In this example, even though the exponential matrix solution stepsize in Fig. 
8 is approximately seven times the central difference stepsize in Fig. 9, the exponential matrix 
solution is much more accurate in predicting the true response. 
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Fig. 8.- Comparison of Exponential Matrix Solution Steps with the True Solution. 



Time (Sec.) 

Fig. 9.- Comparison of Central Difference Solution Steps with the True Solution. 
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Example with Three Degrees of Freedom 

Analogous to the two-degree-of freedom example, consider a cantilever beam with the same 
0.1 m square cross section, 3 meters long, with the three condensed coordinates as shown in Fig. 
10, using the same material constants as in the previous example, the equations of motion for the 
shown coordinates are written as: 
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As in the previous example, first a single step exponential matrix solution is compared with 
the true solution for the free vibration response. The results are depicted in Fig. 11. It is 
seen that increasing the number of degrees of freedom does not diminish the usefulness of the 
exponential matrix solution. However, as may be expected, a more accurate computation of the 
matrix exponential is required to keep the same time increment as the number of degrees of 
freedom is increased. 


DISCUSSION AND CONCLUSION 

The presented simple examples represent only the specific case of undamped free vibrations 
of a structure. However, the demonstrated benefits are also applicable to the analysis of dynamic 
systems under transient loading. In that case using the exponential matrix solution, the solution 
stepsize would depend only upon the discretization requirements of the applied transient loading. 
One could always include a sufficient number of terms in the computation of the exponential 
matrix to satisfy the time-history solution requirements of the high frequency structural response, 
using a time increment of any size. Eq. (16) can be used in the direct integration of the equations 
of motion to increase the practical stepsize and to improve the solution accuracy. Using Eq. (16) 
it is possible to keep the numerical integration stepsize at a practical time increment without 
losing the effects of high frequency structural response. 

The exponential matrix solution may also prove useful for the time-history dynamic analysis 
of nonlinear structures under transient loading. In particular, the exponential matrix method 
would be useful when nonlinear structureal behavior does not significantly affect the high fre- 
quency response. In this case the exponential matrix solution stepsize can be selected to accommo- 
date the nonlinear response characteristics of the structure as well as the accurate representation 
of the transient loading. 

There has been a significant number of publications by computer scientists on the effective 
computation of the matrix exponential [2, 3, 4, 5]. Algorithms using the Pade approximations 
appear to be the most successful ones from a survey of the literature [6,7,8,9,10]. It is possible to 
make a realistic estimate of the accuracy of the fundamental solution at a given time t. Again, it 
is easier to achieve higher levels of accuracy near the origin t = 0. 

The exponential matrix solution should be further evaluated as a tool to improve the solution 
accuracy at a practical stepsize in the direct integration of equations of motion in structural 
dynamics. The stepsize is limited only by structural nonlinearities and the computational accuracy 
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Fig. 10.— Example Beam with Three Coordinates. 
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Fig. 11- Comparison of Exponential Matrix Solution Steps with the True Solution 
for Example Beam with Three Degrees of Freedom. 

(Exponential Matrix Solution is Based on 30 Terms.) 
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of the matrix exponential. The relationship between the accuracy of the matrix exponential and 
the maximum stepsize needs to be quantified in general. 

Having a method that can take larger intervals in the time domain would pave the way for 
more efficient finite time element algorithms to become practical in the simulation of dynamic 
response, enabling the use of different time intervals and/ or different time quadrature rules to 
be used as needed at different locations of the structure and at different points on the time 
coordinate. 
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